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Abstract 

We conduct a numerical study of the dynamic behavior of a dense hard sphere 
fluid by deriving and integrating a set of Langevin equations. The statics of the system 
is described by a free energy functional of the Ramakrishnan-Yussouff form. We find that 
the system exhibits glassy behavior as evidenced through stretched exponential decay and 
two-stage relaxation of the density correlation function. The characteristic times grow 
with increasing density according to the Vogel-Fulcher law. The wavenumber dependence 
of the kinetics is extensively explored. The connection of our results with experiment, 
mode coupling theory, and molecular dynamics results is discussed. 

1992 PACS numbers: 64.70Pf, 61.20Ja, 61.20Lc 
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I. INTRODUCTION 

Despite extensive experimental, numerical and theoretical investigations 1 ' 2 over 
several decades, the present understanding of the slow non-exponential dynamics of dense 
liquids near the glass transition remains incomplete. When a liquid is cooled rapidly enough 
to temperatures below the equilibrium freezing temperature, crystallization is bypassed 
and the system undergoes a transition into an amorphous solid state called a glass. The 
characteristic relaxation time t* , as reflected in a large number of experimentally measured 
quantities such as viscosity and dielectric relaxation, grows rapidly in the supercooled state 
as the temperature is decreased or the density increased. The glass transition temperature 
(or alternatively, density) T g (p g ) is conventionally defined as the temperature (density) 
where the viscosity reaches a value of 10 13 poise. 

In recent years, considerable progress in the development of a theoretical and ex- 
perimental understanding of the phenomena associated with glass formation and the glass 
transition has been achieved. In particular, the so-called mode coupling (MC) theories of 
the glass transition 3-10 have led to a framework for understanding and interpreting 
many experimental results. In MC theories, the slowing down of the dynamics near the 
glass transition is attributed to a nonlinear feedback mechanism arising from correlations 
of density fluctuations in the liquid. The MC equations for the dynamics near the glass 
transition were originally derived 4 ' 5 from the kinetic theory of dense fluids and were later 
generalized and extended by Gotze and co-workers 6,7 . In the original version of MC 
theories 4 ' 5 , the characteristic time scales of the liquid are predicted to exhibit a power 
law divergence at an 'ideal glass transition' or crossover temperature T c , higher than T g . 
However, this divergence is not found experimentally: the power law form breaks down, 
and relaxation times at T c are typically of order 10~ 8 s. More recent calculations 9 ' 10 have 
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uncovered a cutoff mechanism which is supposed to round off the predicted divergence and 
to restore ergodicity over a much longer time scale. Most of the existing MC studies of 
the glass transition do not take into account the equilibrium short-range structure of the 
dense liquid. As a result, these theories do not lead to any prediction for the wavenumber 
dependence of the relaxation. Some attempts towards the incorporation of information 
about liquid structure in the MC formalism have been made recently 8 ' 10 . MC theories 3,7 
and experimental information can be combined in a consistent description of the decay of 
«S(q, £), the spatial Fourier transform of the dynamic density correlation function, at suf- 
ficiently low temperatures or high densities, when the decay times have grown very large. 
The decay according to this scheme takes place in a succession of several regimes: after a 
fast decay in times of order of the inverse phonon frequency, a first slow decay occurs which 
MC theories predict to be an inverse power law in time. This is called the /3-relaxation 
regime. Evidence for power-law decay of correlations in the (3 regime is provided by light 
11 and neutron scattering experiments 12 .This decay is to a nonzero value (an apparent 
nonergodic phase) from which the system eventually moves away leading to the primary 
or a-relaxation regime. The relaxation in the a regime is found to follow the so-called 
Kohlrausch - Williams - Watts 'stretched exponential' form 13 . Both the duration of the 
(3 relaxation and the time scale of the stretched exponential decay in the a regime are 
found to increase sharply as the 'glass transition' is approached. In some cases, the (3 
and the a regimes are separated by a region of so-called von Schweidler relaxation which 
has a power law form. The stretched exponential behavior in the a regime and the von 
Schweidler relaxation are seen in dielectric measurements 14 and in light 15 and neutron 
scattering 16 experiments. The parameters which describe the decay of S(q,t) are known 
to be (/-dependent, but the nature of this dependence has not been studied in detail. Thus, 

4 



the MC theories provide a qualitative understanding of a number of experimentally ob- 
served features of glassy relaxation. However, some of the detailed MC predictions are 
not in agreement with experiments 17 and the MC description clearly fails to account 
for the behavior observed at temperatures close to and lower than T c . It is generally be- 
lieved 3 that this failure arises from the fact that MC theories do not take into account 
activated processes involving transitions between different local minima of the free energy 
which are supposed to develop as the temperature is lowered below the equilibrium freezing 
temperature. 

In contrast to MC theories which portray the glass transition as being purely dy- 
namic in nature, there have been a number of attempts 18 ~ 21 to develop a 'thermody- 
namic' theory in which some of the interesting behavior observed near the glass transition 
is attributed to an underlying continuous phase transition. These attempts have been 
motivated by the fact that the observed growth of the relaxation time in so-called 'fragile' 
liquids 2 is well-described by the Vogel-Fulcher law 22 : 

t*=T 8 e T =*S (1.1) 

where r s is a characteristic microscopic time, a is a dimensionless constant and T is a 
characteristic temperature which is found to be well below the conventional glass transi- 
tion temperature T g . This law suggests the possibility of a so-called 'thermodynamic glass 
transition' that would take place at the temperature To if thermodynamic equilibrium 23 
could be maintained down to this temperature. Such a transition is also suggested by the 
fact 2 that the Kauzmann temperature 24 at which the entropy difference between the 
supercooled liquid and the crystalline solid extrapolates to zero is close to T . Since in 
practice the system falls out of equilibrium at temperatures close to T g , a direct experi- 
mental test of the existence of such a transition is not possible. Also, any calculation that 
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explicitly demonstrates the existence of such a transition in a physically realistic system is 
not yet available. Thus, the 'thermodynamic glass transition' scenario remains essentially 
speculative. In this approach, the growth of the relaxation time is attributed to a growing 
correlation length which would diverge at the ideal glass transition. Two recent numerical 
studies 25 ' 26 which looked for such a growing correlation length did not find any evidence 
for its existence. A recent experiment 27 , on the other hand, has presented evidence for 
the existence of a length scale that grows as the temperature is lowered below the crossover 
temperature T c . 

The static and dynamic properties of dense liquids have also been studied exten- 
sively by molecular dynamics (MD) simulations 28 ~ 32 . The very nature of the simulation 
method restricts such studies to simple model systems and to time scales which are rather 
short (of the order of 10~ 9 sec). In spite of these limitations, MD simulations have pro- 
duced a number of interesting results which are in qualitative agreement with predictions 
of MC theories and results of experiments (such as those using the neutron spin-echo tech- 
nique 33 ) which have time scales comparable to those of the simulations. This observation 
leads to the interesting and useful conclusion that a study of simple model systems may 
be sufficient for understanding the basic physics of the glass transition. 

It is evident from this brief survey of the current status of the glass transition 
problem that an obvious need exists for the development of new analytic and numerical 
methods which may address some of the outstanding issues related to this problem. In this 
paper, we describe the results obtained from the application of a new numerical method to 
a study of the dynamic behavior of a dense hard-sphere liquid near the glass transition. The 
method we use consists of direct numerical integration of a set of Langevin equations which 
describe the nonlinear fluctuating hydrodynamics (NFH) 34 of the system. Information 
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about the static structure of the liquid is incorporated in the Langevin equations through a 
free-energy functional which has a form suggested by Ramakrishnan and Yussouff (RY). 35 
Recently, it has been shown 36 that the RY free energy functional provides a correct 
mean-field description of the statics of the glass transition in this system. In that work, 
a numerical procedure was used to locate local minima of a discretized version of the 
RY free energy appropriate for the hard sphere system. A large number of glassy local 
minima with inhomogeneous but aperiodic density distribution were found to appear as 
the average density was increased above the value at which equilibrium crystallization 
takes place. At higher densities, the free energies of these minima were found to drop 
below that of the minimum representing the uniform liquid signaling a mean-field glass 
transition. The success of the RY free energy functional in providing a correct description 
of the statics of the glass transition of the hard sphere system suggests that a good starting 
point for a study of the dynamics of this system would be obtained by incorporating this 
free energy in the appropriate NFH equations. 

A number of important issues are addressed in our study of the dynamics. By 
comparing the results of our calculations with existing MD results 28 on the same system, 
we are able to test the validity of the NFH description which is cast in terms of coarse- 
grained number and current density variables instead of the coordinates and momenta 
of individual particles. The correctness of the NFH equations we use, although usually 
taken for granted, is not obvious in view of the fact that the hydrodynamic terms in these 
equations describe the physics at relatively long length scales, whereas the terms arising 
from the free energy functional involve length scales of the order of (or smaller than) the 
interparticle spacing. In the RY free energy functional, information about the microscopic 
interactions is incorporated in the form of the Ornstein-Zernike direct pair correlation 
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function 37 of the liquid. This appears to be adequate for a correct description of the statics 
of the freezing of the liquid into both crystalline 35 and glassy 36 states. One of the questions 
we address in the present study is whether this is also sufficient for a correct description of 
the dynamic behavior. Two numerical studies of NFH equations describing the dynamics 
of dense liquids have been reported recently 38 ' 39 . The main difference between these 
studies and the present one is that the equilibrium structure of the liquid was treated 
only in an approximate way in the earlier calculations. Consequently the results of these 
calculations could not be compared directly with MD data. The results of these calculations 
did not exhibit some of the features (such as two-regime decay of correlations) generally 
associated with glassy dynamics, although non-exponential decay of correlations often 
associated with its onset was found. In Ref. (39), it was suggested that this failure arises 
from the approximations made in the treatment of liquid structure. The present study 
provides the opportunity to check whether this explanation is correct and to investigate in 
general the role of static structure in the dynamics of a dense liquid. Further, perturbative 
treatments of the NFH equations similar to the ones we use are known 9 ' 10 to lead to results 
which are very close to those obtained from the MC approach. Apart from numerical 
errors arising from spatial discretization and the integration procedure, our treatment 
of these equations is exact. In particular, our numerical solution of these equations is 
obviously nonperturbative, therefore, a comparison of the results of our calculation with 
MC predictions provides a way to test the validity of some of the approximations made in 
the analytic studies. Since our calculation takes full account of the short range structure 
of the liquid, it has the most direct connection with the MC studies of Kirkpatrick 8 and 
Das 10 . Finally, by monitoring which minima of the free energy are visited during the time 
evolution of the system we are able to determine whether the observed dynamic behavior 
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arises from nonlinearities of density fluctuations in the liquid or from transitions among 
different glassy minima of the free energy. It is not possible to distinguish between the 
effects of these two kinds of processes in conventional MD simulations. 

We have studied the time-decay of S(q,t), the spherically averaged spatial Fourier 
transform of the time-dependent density correlation function of the hard-sphere liquid, 
for different values of wavevector q and for a number of values of the reduced density n* 
(n* = pod 3 , where po is the average number density and a is the hard sphere diameter) 
in the range 0.75 — 0.93. It was found in Ref. (36) that a hard sphere system described 
by a discretized version of the RY free energy exhibits a crystallization transition near 
n* = 0.83. Since the present calculation uses the same discretized free energy as that of 
Ref. (36), the system may be considered to be in the 'supercooled' regime for a large part 
of the density range considered by us. The numerical efficiency of the Langevin dynamics 
arising from the use of coarse grained variables enables us to verify that the statics of the 
system remain stationary throughout the long time intervals that we consider. The ob- 
served dynamic behavior described later in detail exhibits a number of characteristic glassy 
features. These include stretched exponential decay of correlations, two-stage relaxation, 
and Vogel-Fulcher growth of relaxation times. Our results are in agreement with existing 
MD data 28 on the dynamics of the hard-sphere liquid and in qualitative agreement with 
other MD results obtained for similar but different systems. This observation establishes 
the correctness of the NFH description used in this work and also demonstrates that the 
RY free energy contains the essential physics of the dynamics of this system. Our calcu- 
lations also reproduce qualitatively a number of predictions of MC theories, however, we 
find significant deviations from some of the quantitative MC predictions for the glassy 
kinetics of the hard-sphere system. Our study indicates that the onset of glassy features in 

9 



the decay of S(q, t) occurs at relatively lower densities for wavevectors close to the first and 
second peaks in the static structure factor. This result clearly illustrates the important role 
played by the equilibrium structure in the long-time dynamics of the liquid. The observed 
g-dependence of the decay of S(q,t) also suggests that the glassy behavior sets in earlier 
(at lower densities) at shorter length scales. Finally, the system is found to fluctuate about 
the liquid-state minimum of the mean field free energy in all our simulations. 

The rest of the paper is organized as follows. Section II contains a description of 
the model considered by us. We first define the model, discuss its statics and derive the 
appropriate NFH equations. By defining appropriate units of length, time and mass, these 
equations are then written in dimensionless form. The method used by us to integrate 
these equations forward in time is described in Section III. We also discuss in this Section 
the physical quantities measured, the data collection procedure, tests of equilibration and 
other related matters. The results obtained from the numerical work are described in 
detail in Section IV. We also compare and contrast our results with those obtained from 
MD simulations of the hard-sphere and similar systems and the predictions of MC theories. 



We have explained in the Introduction that we will analyze the numerical solution 
of a set of Langevin equations appropriate to a dense hard-sphere fluid. We begin here by 
discussing the statics of the model. These are given in terms of a free energy which is a 
functional of the two fields in the problem: the density field p(r, t) and the current field 
g(r, t). This free energy has two terms: 



II. THE MODEL 




(2.1) 
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where F[p] is of the RY form 35 : 

F[p] = Ftipo] + k B T (^J dr(p\og(p/p ) - dp) - (1/2) J dvdv'C{v - r') (2.2) 

In (2.1) and (2.2) p is the number density field and 5p = p — po the deviation of that field 
from its average value po- Fi is the free energy of the uniform liquid, T is the temperature 
and ks the Boltzmann constant. In the first equation mo denotes the mass of a hard 
sphere. Finally, C(r — r') is the direct correlation function. The inclusion of this function 
in the free energy ensures that upon linearization of the logarithm in the first term on the 
right hand side of (2.2) one obtains the usual expression for the static structure factor of 
a simple fluid in terms of C. For hard spheres a simple expression for C(r — r') can be 
obtained in the Percus-Yevick approximation 37 : 

C(0 = -Ai - QrifXit - (1/2^/Aie 2 ; £<1 (2.3a) 
C(0 = ; C>1 (2-36) 

where £ = |r — r'\/a, rjf is the packing fraction: 

r?/ = (7r/6)p ^ 3 = (tt/6) n* (2.4) 

and: 

Ai = (1 + 2?7/) 2 / (1 - rjf) 4 (2.5) 

A 2 = -(l+r ?/ /2) 2 /(l-r ?/ ) 4 (2-6) 

We have written po in the denominator of the g 2 term in Eq.(2.1), rather than p as in 
Ref. (39), so that the full density dependence of the free energy is given by the RY form. 
As pointed out in Ref. (41), if we used p we would obtain upon functional integration 
with respect to the gaussian variable g 2 a log(p) contribution involving density fluctuations 
associated with the kinetic energy. These are already included in Eq. (2.2). 
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It is convenient to choose units at this point, which will simplify subsequent expres- 
sions. We take mo as the unit of mass. We choose also a unit of length, h, which we shall 
later identify with the lattice constant of the computational lattice. For the unit of time 
we choose to such that: 

t = h/c (2.7) 

where c is the speed of sound. This choice of the unit of time is motivated as follows: the 
characteristic phonon time for the system is t p = l/(cq) where q is a wavevector. Hence 
tp/to = l/(qh). We shall choose below the ratio a/h so that qh is of order unity in the 
wavevector range of interest. Hence the choice (2.7) corresponds to a characteristic time 
of order t p . The same choice was taken in Ref. (39). 
We then define the dimensionless quantities: 

x = r/h (2.8) 

n = ph 3 (2.9) 
j = gh 3 /c (2.10) 
and introduce also the dimensionless free energy F[n, j] in terms of the above quantities: 

F[n,j] = (1/2) / d* 3 — + F n [n] (2.11) 
J no 



F n [n] = F ln [no] + K (^J dx(n log(n/n ) - Sn) - (1/2) J rfxrfx'C(x - x.')SnSn^j (2.12) 



where 



K=^L (2.13) 
moc z 
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The quantity K is for hard spheres a function of only the density. This follows from: 

c 2 = (2.14) 
m po« 

where the compressibility k is related to the structure factor S(q) by ksTpoK = S(q = 0). 
The latter quantity can straightforwardly be calculated from (2.3) and one obtains: 

K=(l + (47r/3)n*(A 1 + (9/2) Vf X 2 + ( Vf /4)X 1 )y 1 (2.15) 

which indeed depends on the density only. The quantities Ai, \2,Vfi an d n* were defined 
above. 

We quote here the relation between our unit of time to and the Enskog mean collision 
time t c commonly used in MD work 42 on hard spheres. One easily finds that: 

t- 1 = 4y^n*g(a)\^(h/a)to 1 (2.16) 

where g(a) is the pair correlation function at contact 43 . The static properties of a 
system described by the RY free energy, and its relation to glassy behavior have been 
discussed in the literature 36 ' 35 and in the Introduction. We turn now to the dynamics. 
The hydrodynamic variables are of course the fields n and j. We denote, following the 
ususal notation, this set of Langevin variables as {ip a } = {^(x, £), j(x, £)}. Formally, the 
Langevin equations can be written as 9 ' 44 : 

where V is the matrix of transport coefficients, 6 are the gaussian noise fields, and the V a 
are the streaming velocities, which incorporate the nondissipative part of the equations of 
motion. Proceeding as in Refs. (9) and (44) one obtains: 

K = -£v>(x)— ] (2.18) 
13 



V * = - E - E*(W^ (2.19) 

J J 

where F and F n are defined in Eqns. (2.11) and (2.12). Following then on precisely as in 
Ref. (39) we obtain, after straightforward algebra the equations of motion: 

<9n(x, t) 



dt 

and: 



+ (l/no)V(ry) = (2.20) 



Ji = ~nV~ - (l/n )J2^jUdj) ~ (1/no) ^jjVijj + (l/n )v^ 2 ji + (2-21) 



dt 5n 

3 3 



where n is the bare shear viscosity in our units 45 . The noise fields 0j(x, t) satisfy the 
second fluctuation-dissipation theorem in the form: 

< e i (x,t)0 J (x / ,t / ) >= -2KXr]n S l}J V 2 S(r - r')5(t - t'). (2.22) 



We recall that in these equations, and in the remainder of the paper, the time is 
measured in units as given in (2.16). In (2.22) the angular brackets denote the thermody- 
namic average. The quantity A is a dimensionless measure of the equilibrium fluctuations. 
The average value uq is related to n* = na 3 through no = n*(h/a) 3 . For hard spheres, 46 
one can write rj in terms of K and the density in the form: 

T} = (V^) 2 v / ^[5/(16 v / ^)][^(a)- 1 + 0.8(27m*/3) + 0.761(27m73) 2 #(a)] (2.23) 

The equations of motion (2.20) and (2.21) are slightly more complicated than the cor- 
responding equations for the model in Ref. (39). There is an additional convection term 
in the second equation and some additional factors of n/n . These complications can be 
directly traced down to the different density dependence of the first term in (2.11) as dis- 
cussed above, that is, to the fact that a kinetic energy contribution is now included in the 
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first term of F n . These differences then arise because of the different form of 8F/5j which 
results in more complicated expressions for the streaming velocities. Our expressions cor- 
rectly reduce to the continuity and Navier-Stokes equations in the linear limit. In general, 
one can interpret (2.20) as the continuity equation if the field g is thought of as po times 
the velocity. The current density is then not given simply by g, however. These questions 
do not affect, obviously, the validity of the conclusions obtained from this model. 

A salient feature of Eq.(2.21) of considerable importance is that the term in 5F n /Sn 
which involves the direct correlation function C is an integral over space with range a. 
Thus, we have to solve in this case not merely a set of partial differential equations with 
stochastic terms but as far as the spatial dependence is concerned an integrodifferential 
equation. This complication makes this a difficult problem to solve numerically. The 
methods used will be explained below. 

III. METHODS 

In this Section we discuss the methods we use to solve our equations of motion, the 
physical quantities we focus on, the data collection methods and statistics, the range of 
parameter values we have studied, and related matters. 

Our objective is to study the dynamic correlations of our system as defined below. 
In order to do so we solve numerically Eqns. (2.20) and (2.21) on a three dimensional 
cubic lattice of size N s . The main complication we must consider is the spatial integral 
involving the direct correlation function C(r). For each integration step this integral must 
be done N 3 times since it is computed over a sphere with its origin at each one of the 
lattice sites. To avoid repeated calculations we create, in the initialization of our program, 
a table listing for each lattice site the location of all neighboring sites to be integrated over 
and their corresponding values of C(r). Additionally, since the sphere is imbedded on a 
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coarse discrete lattice we use a finer mesh than defined on the original lattice to improve 
the accuracy of the integration. The procedures employed to integrate the remaining set 
of differential equations over time and generate the gaussian noise are identical to those 
used in Ref. (39) and references cited there. 

We will focus our analysis on the time dependence of the dynamic structure factor 
S(q,f): 

5(q, t)= j dW q ' (x - x ' } < 6n(x, 0)<5n(x', t) > (3.1) 

Specifically, we will consider the angular average of S(q,t). On a cubic lattice, it is 
appropriate to define the effective length of q, q, as: 

q 2 = 2(3 — cosq x — cosq y — cosq z ), (3.2) 

and we perform angular averages of q-dependent quantities by averaging over values of q 
in the first Brillouin zone, in a spherical shell of mean radius (as given by q) corresponding 
to that of the vector (ttQ/N, 0, 0) and thickness w/N. The value of Q ranges from 1 to 
approximately 3 1 ' 2 N, although only a smaller range is free of finite size effects. We will, 
for simplicity of notation, denote quantities averaged in this way by simply dropping the 
vector symbol from the wavevector argument: S(q, t), and often we will indicate the values 
of q by the 'shell number' Q. 

Next we turn to our choices for parameter values. The correlation functions that we 
are interested in are spatially short ranged. It is therefore not necessary to use extremely 
large lattice sizes. The results presented were obtained using N = 15. As in Ref. (39), 
this proved to be adequate. We checked that finite size effects do not affect the dynamics 
in the wavevector region for which results are presented here, (5 < Q < 15, see Section 
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IV) by performing a portion of the calculations (with reduced statistics) at N = 25. Our 
choice of the ratio a/h which fixes the length scale for the problem, is dictated by two 
concerns. The first is that we wish to be able to study the dependence of the dynamics 
on wavevector in the main region of interest from the point of view of the static structure 
factor S(q) = S(q,t = 0). Thus, we wish to choose our unit of length so that the main 
peak in S(q) falls in the middle part of the range of wavevectors within the first Brillouin 
zone of the computational lattice. Secondly, to avoid crystallization at the higher densities 
studied, we have found that we need N and a to be incommensurate. Selecting a/h = 4.6 
leaves q m ax well away from the zone edge for all densities, near the Q = 8 shell, and is 
clearly incommensurate with N. The density range we have investigated includes n* = 0.5, 
0.75 < n* < 0.90 at 0.05 intervals, and n* = 0.93. The first of these is well within the 
dilute liquid region, and was used only to check that limit. As explained in Ref. (39), it is 
necessary to include the parameter A in order to represent the actual fluctuations through 
gaussian noise. Its precise value is not crucial, since it essentially amounts to a choice of 
the normalization of the static fluctuations S(q), but it clearly must be small since the 
density must always be positive. We have taken here A = 0.001. The remaining parameters 
K and r\ are functions of n* and may be calculated as discussed in the previous Section. 

We turn now to the very important question of data collection. We are particularly 
concerned with ensuring that within statistical error the averages collected be equilibrated, 
that is, stationary in the time scales studied. We find that, with the data collection 
procedure as outlined below, the initial conditions that we use to begin the integration of 
the equations of the motion are unimportant (except in that they determine to some extent 
the duration of transient behavior) and we usually take them to be a flat distribution of n 
equal to its average value, no, and vanishing currents. We then monitor the current-current 
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correlations as a function of running time to- After a relatively short time of order 10 the 
current correlations reach their equilibrium value as given by the equipartition theorem. It 
might be tempting to assume that the density correlations have also equilibrated by that 
time and such an assumption is sometimes made in MD work, but we find that for our 
system at least this assumption does not hold. 

To study the density correlations, we store, for running times to > tx, where to is 
the time measured from the initiation of the computation, at a large number of periodic 
time bins, the products of the form <5n(x, to)Sn(x.', t + i) for all x x'. We then monitor 
the spherically averaged spatial Fourier transform, S(q,t,to) of the quantity: 

5(x, x', t, t ) =< Sn(x, t )6n(x!, t + t)> (3.3) 

where the average is understood to be over a number n& of time bins separated by an 
interval At. The time range covered by the averaging process is tn = rib At. In order 
for S(q,t,to) to be an adequate approximation to the thermodynamic average S(q,t) it is 
required that it be not only independent of to, but also independent of tn within statistical 
error. Dependence on the to indicates the presence of a transient. Dependence on tn 
indicates that the averaging time is too short for ergodicity to hold. A very important 
point is that we find that the minimum value of the transient time for density fluctuations 
is not t K , but it is of the order of the slowest characteristic decay time t* in the system. As 
discussed in the next Section, t* is a strongly increasing function of density, and is much 
longer than the equilibration time for the kinetic energy. Similarly, it is necessary for the 
average to include a range tn of order of several times t*. We find that S(q, t, to) at higher 
densities has considerable oscillations over t time ranges smaller than t*. As one increases 
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the density an estimate for t* can be found by extrapolation from the lower densities, since 
the behavior of t* with density turns out to obey the Vogel-Fulcher law 22 . 

Working within these constraints, obtaining statistically reliable results still requires 
averaging over a large number of time bins which means very large total running times. 
In addition, in order to eliminate any possibility of spurious correlations due to a peculiar 
transient, we have repeated the whole procedure three to five times at each density. The 
results presented here correspond to a combined total of between 1000 and 3900 bins, 
depending on the density, at all densities we present results for except n* = 0.75 where 
a total of 600 was taken. These very large numbers, much larger than the corresponding 
numbers in Ref. (39), should be considered as comparable to the 'number of runs' in a 
standard simulation for a nonequilibrium problem such as spinodal decomposition and 
lead to very good quality data. The cost of obtaining the data rises accordingly, of course: 
a total of 200 hours of Cray2 and Cray X/MP time were required. 

We have also verified that the quantity S(q, t = 0) calculated following the above 
procedure and over the time ranges just described is consistent with the purely static 
result. To do this, first we calculate the static result: we evaluate numerically, using fast 
Fourier transforms, the discrete Fourier transform C(q) of (2.3) for a system of the size 
considered, and we obtain from that result the static result S s (q) oc 1/(1 — n*C{q)) . To 
obtain this equation one must expand the first term in (2.12) to second order in 5n. The 
comparison between S(q) and S s (q) is shown in Fig. 1. We can see that the two results 
are in very good agreement at the densities plotted. The PY static values are known to 
be (see e.g. Fig 2 in Ref. (42)) in good agreement with MD results in this density range. 
Our results are then also in agreement with MD in this limit. At higher densities, the 
computed result represents a higher degree of order than that calculated from the statics. 
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This is due to the fact that the expansion of the free energy just alluded to is not as well 
justified at those densities. When averaged over time scales of less than several times t*, 
S(q : 0) oscillates broadly about its stationary value. Our results for S(q, 0) at all densities 
studied show that we are dealing with a liquid-like state here. If we use a commensurate 
value of cr, or if we increase n* to 0.95, we find indications that crystallization begins. We 
plan to study this crystallization question in future work. 

The stability of our dynamically obtained results over long computational runs, and 
their agreement with static results constitute a very stringent check of the stability of our 
numerical algorithms. 

To analyze our results, it is convenient to introduce the normalized quantity C(g, t) 
defined as: 

where, we recall, we are dealing with angularly averaged quantities as explained above 
with q defined in (3.2). We will then, in the next Section, characterize the decay for this 
quantity for all values of q. By averaging C(q, t) over a large effective number of runs we 
are able to obtain results sufficiently smooth to confidently fit functional forms to the data 
as the next Section demonstrates. 



IV. RESULTS AND DISCUSSION 

We now turn to the discussion of the normalized dynamic structure factor (see 
Eq. (3.4)), C{q,t). We have first verified that at low densities (n* = 0.5), this quantity 
decays exponentially at all wavevectors studied. We will consider in what follows only 
the more interesting range 0.75 < n* < 0.93. We include in all of our fitting attempts 
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all data up to either the maximum time for which we have data at the particular values 
of n* and q under consideration or up to the time where C(q, t) is so small that it fades 
into the noise. This occurs when C{q,i) < 0.025 except in some very few cases where the 
data becomes noisy at the 0.05 level. As in Ref. (39) the statistical noise seems to have 
an additive component which causes the relative errors to increase when C(q,t) becomes 
smaller. We recall that the relation between our time units and Enskog collision time is 
given by (2.16). The factor relating the two inverse times varies from 1.52 at n* = 0.75 to 
1.90 at n* = 0.93. The relation between our unit of length and a is the trivial factor of 
a/h = 4.6. 

We begin by attempting a fit to a stretched exponential form: 

C(q,t) = e - (t/To) " (4.1) 

where the parameters tq and (3 are functions of n* and q. This step is motivated in part 
by the expectation from preliminary inspection of the data that there is a wide region of q 
and n* values for which this form is adequate. Indeed we find that for a considerable part 
of the data particularly in the lower density region this turns out to be a satisfactory fit. 
Even in the cases where the data is not well fitted by the form of Eq. (4.1), the number 
To (q, n*) still gives a useful figure of merit or overall estimate of the decay time. The results 
of this fit are in Tables I and II. We have indicated in these Tables the cases where the fit 
to the above form is a good fit to the data and those in which it is actually not the best fit 
by enclosing the latter cases in parentheses. The quality of the fits can be judged visually 
and by the x 2 values that we obtain. 

The salient points of the results in Table I are apparent: the characteristic time is 
at constant density a very strong function of q, and it has a maximum at the wavevector 
shells close to where the static structure factor S(q) has its maximum. A less well defined 
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secondary maximum at wavevectors close to the second maximum of S(q) can also be 
discerned. The characteristic time is not however of the simple form To(q, n*) oc i]S(q) as 
one might guess naively: it is neither proportional to S(q) at constant density (recall r\ 
depends on the density only) nor proportional to r\ at constant q. This indicates as might 
be expected, a q-dependent renormalization of the transport characteristic times. This is 
more marked at higher densities. The stretching parameter f3 (Table II) remains relatively 
close to unity, although stretching is evident as one increases the density particularly at 
shorter wavelengths. We note that the entries in the Table corresponding to the smaller 
values of (3 are in most cases not good fits to the data, as we shall discuss below. 

It is apparent from our tabulated results that the slowest overall decay is found, in 
our spherically averaged results, at the shell number Q = 8, near the main peak in S(q). 
We consider then the main decay time To at that wavevector, as representing the slowest 
transport time at that density t*. We plot in Fig. (2) this quantity as extracted from Table 
I as a function of density. We have fitted this quantity to the Vogel-Fulcher 22 law (1.1) 
which when the density, rather than the temperature, is the externally controlled variable, 
takes the form: 

t*(n*) =ae~ l/{v ~ Vc) (4.2) 

where v = 1/n*. As shown in the Figure, we find an excellent fit with a value of v c 
corresponding to n* = 1.23. This is in excellent agreement with the analysis of MD data 
of Ref.(28) where the result n* = 1.21 was obtained. A three parameter fit to a power law 
is also adequate, but the quality of the fit is not quite as good. 

We now turn to the decay modes for which the stretched exponential is not a 
satisfactory fit. We have examined the data very carefully and attempted a large number 
of fits to different proposed decay forms. We find that the reason for the failure of (4.1) 
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is usually that the decay is a two step process, or in some cases that the data appear to 
decay to a finite constant in the time range studied. Accordingly, we have found that the 
best form that fits this portion of the data is : 

C{q,t) = {l-f)e-W+fe- t " (4.3) 

where it is understood that the time scales are well separated, r' 3> r. One can also include 
a stretching exponent in the second decay, but we have found it usually unnecessary and 
we were leery of multiplying the number of fitting parameters. We note that the form 
(4.1) is a particular case of (4.3) with / = or, alternatively, r ~ r' . It also represents a 
decay to a constant (quasinonergodic behavior) whenever r' = oo, which of course should 
be understood as meaning that the second decay time is beyond the region fitted. Finally, 
with / finite, (4.3) yields a region of apparent power law decay of the von Schweidler 
form when r <C t <C r' and the time scales are widely separated. We find that this 
very general form fits well (again, as measured by the appropriate reduced x 2 values and 
visually) all cases where as indicated by the parentheses in Tables I and II the single 
decay fit was not satisfactory. The parameters r, r', and / are given in Table III. The 
parameter f3 turned out to be very close to unity in all of these cases. In Figures (3) and 
(4) we show a wide selection of our results and the corresponding best fits, according to 
the parameters in Tables I-III. For every density we have two Figures, which display the 
decay of C(q,t) at that density for a sampling of wavevector values. The values shown 
are Q = 6, 8, 10, 12, 13, 14. These include all Q values where the decay is not of the form 
(4.1). Values skipped have been omitted to avoid excessive clutter. Keeping in mind the 
Figures and the three Tables, we see very clear trends. First, as one increases the density, 
nonexponential behavior crops out first at relatively larger wavevectors, that is, shorter 
distances, which is physically sensible and was observed also in Ref. (39). Specifically, this 
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behavior appears first at values beyond the main peak in the structure factor, in the region 
of the second peak. As the density increases further the beginning of the same phenomenon 
is observed at smaller wavevectors at or near the main peak in S(q). This phenomenon 
appears to be related to the q dependence of (3 found in mode coupling calculations 47 . The 
characteristic times increase strongly with density, and so does their relative separation. 
They are also strongly q dependent. On the other hand, / appears to depend rather weakly 
on n* but more strongly on g, and its q dependence appears also to correlate with that 
of the statics. From the point of view of (4.3) it appears that the lower density stretched 
exponential form (4.1) should be viewed as arising from a confluence of the time scales, 
rather than from the vanishing of /. 

We monitor also the value of the mean field free energy during the computation. We 
find that the system fluctuates around a liquid like state minimum. Thus, our simulations 
establish that the first three or four orders of magnitude in the growth of the relaxation 
time in the hard sphere system arise from the effects of nonlinearities associated with the 
realistic free energy used. 

It is very instructive to discuss our results in connection with MD work, with mode 
coupling theoretical work on glassy systems, and with experiment. Some caution is re- 
quired, however, since it must be kept in mind that we are dealing here with relatively low 
densities: the degree of 'supercooling' that our choices of N and a achieve is limited in 
comparison with what can be done experimentally. We have not used other procedures, 
such as mixing spheres of two different sizes, which might allow further increases in the 
density. There are a variety of decay regimes and forms reported in the literature for glassy 
systems. A scenario in favor of which some consensus has gathered 17 is described in the 
Introduction. We believe that our results at higher densities and shorter distances, in the 
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region where the full form (4.3) is obeyed might be interpreted within the above scheme 
as representing a first decay which combines phonon and j3 -relaxation, leading eventually 
to the second decay which might then be identified with a -relaxation. In some cases only 
the quasinonergodic regime is reached in our computations. The von Schweidler regime 
can arise, as explained above, as a precursor of the second decay. The phonon assisted 
decay does not appear as a separate regime in our data, which is possibly due to the coarse 
graining nature of Langevin dynamics with respect to microscopic times. We have verified 
that an inverse power law does not fit the earlier portion of our data as well as a stretched 
exponential form. 

We have already mentioned the agreement of our results for t* and S(q) with the 
MD results of Refs. (28) and (42) respectively. Comparison with MD work on glass 
dynamics is complicated by the fact that such MD results are for different systems, such 
as binary soft sphere mixtures (see e.g. Refs. (30) and (32)) or Coulombic systems (31) . 
Obviously, detailed comparisons will have to await further work. However, as results for 
different MD systems are qualitatively consistent with each other, a qualitative comparison 
with our results is possible. Although comparison of MD and Langevin time scales is in 
general a vexing problem, for our purposes here it might be adequate to simply assign 
to our pho non-based time unit a value of order 0.1 ps. We then find that the results of 
Ref.(30) correspond to a shorter time range than our higher density results, while those 
of Refs. (32) and (31) are comparable. Neither Ref.(32) which covers higher densities than 
ours, nor Ref.(31) follow C(q,t) until it decays to zero or to a small value, as we do, 
and they assume that their systems have reached equilibrium when the kinetic energy 
equilibrates, an assumption which we checked and we found, for our dynamics, erroneus. 
We do not know if this is also the case for MD. Despite all of these caveats, our results 
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are consistent with those of both Refs. (32) and (31). They observe first a glitch at very- 
early times, which they attribute to phonons, and then two relaxation modes separated 
by a slow region, which they identify with (3 and a- relaxations. They do not attempt to 
quantitatively characterize the (3 regime by either a power law or an exponential, but they 
do obtain an exponential form (stretched at higher densities) for the second decay. As 
stated above, our data does not show any early time glitch, which could be due either to 
Langevin coarse graining, or to the glitch being an artifact of poor system equilibration. 
Except for this, our results are consistent with MD in that a naive extrapolation of our 
results to higher densities would predict that the time scales r and r' in (4.3) become more 
separated, and possibly larger values of / appear. We have verified that the MD data of 
Ref.(31) can in fact be fit to the form (4.3) with appropriate parameter values. Thus we 
conclude that our results are consistent with MC as well as MD. 
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FIGURE CAPTIONS 

Fig. 1. The static S S (Q) calculated as explained in the text, compared with the numerical 
results for S(Q) = S(Q,0). Both quantities are spherical averages plotted vs. wavevector 
as given by the shell number Q defined in the text. The curves shown are for n* = 0.75 
(S S (Q) solid line, S(Q,0) medium dashes) and for n* = 0.80 (static results short dashes, 
dynamic results shown as dots). 

Fig. 2. The characteristic decay time, t* plotted as a function of density n*. Data is fitted 
to a Volger-Fulcher form as explained in the text. 

Fig. 3. The normalized correlation function C(q, t) plotted vs. t for Q= 14 (top curve), 12, 
and 6 (bottom curve) . The data and the best functional fits (smooth curves) as explained 
in the text are shown for (a) n* = 0.75 , (b) n* = 0.80 , (c) n* = 0.85, (d) n* = 0.90 , and 
(e) n* = 0.93 . 

Fig. 4. As in the previous figure but for Q— 8 (top curve) ,13, and 10 (bottom curve) at 
t = 100. 
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TABLE CAPTIONS 

Table I.- The parameter To of Eq.(4.1) as a function of density as given by n* and wavevec- 
tor as given by shell number Q. The parentheses indicate data sets for which the fit is not 
satisfactory. These cases are fitted by Eq. (4.3) and Table III. 

Table II.- The parameter f3 of Eq.(4.1) as a function of density as given by n* and 
wavevector as given by shell number Q. 

Table III.- The parameters r,r', and / of Equation (4.3), for the appropriate range of n* 
and Q values, (see text and Tables I and II). The values of r' indicated by infinity must 
be interpreted as being beyond the time window fitted. The three daggered values of r' 
indicate cases where an exponent j3' different from unity was required. In these cases (3 
was set to unity, so the number of parameters was the same, but in the other cases the 
unrestricted (3 remained close to unity. 
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Table I 



0\n* 


n 75 
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f) Q3 
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14 


(134) 


(150) 


(247) 


(221) 


(301) 
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Table II 
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14 
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(0.47) 
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(0.32) 
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Table III 
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